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We test a subgrid-scale spectral model of rotating turbulent flows against direct numerical sim- 
ulations. The particular case of Taylor-Green forcing at large scale is considered, a configuration 
that mimics the flow between two counter rotating disks as often used in the laboratory. We per- 
form computations in the presence of moderate rotation down to Rossby numbers of 0.03, as can 
be encountered in the Earth atmosphere and oceans. We provide several classical measures of the 
degree of anisotropy of the small scales of the flows under study and conclude that an isotropic 
model may suffice at moderate Rossby numbers. The model, developed previously (Baerenzung et 
al., Phys. Rev. E 77, 046303 (2008)), incorporates eddy viscosity that depends dynamically on the 
inertial index of the energy spectrum, as well as eddy noise. We show that the model reproduces 
satisfactorily all large-scale properties of the direct numerical simulations up to Reynolds numbers 
of ~ 10 4 and for long times after the onset of the inverse cascade of energy at low Rossby number. 



I. INTRODUCTION 

Rotating flows are commonplace in nature, the influ- 
ence of rotation being measured by the Rossby number 
Ro = Uq/2LqQ, with Uo the r.m.s. velocity, Lq a char- 
acteristic lengthscale of the flow and the rotation rate. 
The Rossby number of the atmosphere is ~ 0.1 and in 
the ocean it can be as small as 10 -3 . Assuming a con- 
stant rotation rate, the Coriolis force that appears in 
the equations leads to the emergence of wave motions 
which, at small enough Rossby number, can be thought 
as dominating the dynamics. However, at high Reynolds 
number, Re = UqLq/v, with v the viscosity, turbulent 
eddies interact with waves, and inertial waves interact 
nonlincarly (in particular through resonances), so that 
the dynamics become complex. The Rossby number can 
be viewed in this way as the ratio of the characteristic 
time of an inertial wave, tw ~ 1/fl to the characteristic 
time of an eddy, or eddy turn-over time, t^l ~ Lq/Uq; 
when small, the waves are rapid and may dominate the 
dynamics. 

Many studies have been devoted to the exploration 
of rotating turbulence, experimental as well as numer- 
ical and theoretical (see e.g.pj). One expects the flow 
to become quasi bi-dimensional under the influence of 
strong rotation but recent studies show that the dynam- 
ics is more subtle, with three-dimensional eddies possi- 
bly prevailing at small scales. The case of small Rossby 
number can be studied using anisotropic extensions of 
closure models, such as the Eddy Damped Quasi Nor- 
mal Markovian approximation (EDQNM hereafter). In 
such approaches, the closure is obtained by modeling the 
damping of fourth-order cumulants (non-zero for a non 
Gaussian field) by a term linear in third-order moments; 
dimcnsionally, the constant of proportionality is the in- 
verse of a time, or a rate fx, taken in EDQNM to be the 
rates known to be significant in the physics of the prob- 



lem. In the simplest case of non-rotating isotropic and 
homogeneous turbulence, these rates are proportional to 
the inverse of the eddy turn over time ti = l/Ut and of 
the viscous time Trj ~ i 2 jv, expressed in terms of the 
scale £ and the velocity at that scale Ui. In the rotating 
case, the wave frequency becomes relevant as well (see 
3 for an early realization of this concept) , and because of 
the anisotropic dispersion relation of inertial waves, the 
model becomes anisotropic itself in terms of a spectral en- 
ergy distribution that is a function of the wavenumbers 
k±_ and k\\ , where _L and || refer to directions relative to 
the rotation axis. 

Other approaches include weak turbulence theory [|[, 
following the original methodology of Benney and Newell 
[f| , and other resonant wave theories 0, 0] recently shown 
to correspond to an asymptotic limit for flows with wave 
dynamics @. 

The link between resonant theories and the EDQNM 
closure for rotating flows has been analyzed in detail re- 
cently ( 0, [l"o| and references therein) ; it can be simply 
said here that when the global damping rate li (omitting 
triadic scale dependence) is dominated by the inertial fre- 
quency, in the limit of tw —> 0, the fourth-order cumulant 
becomes negligible, the flow becomes quasi-gaussian and 
the closure occurs naturally. Of course, this limit needs 
to be taken carefully, in particular when approaching the 
slow manifold (k\\ = 0) @, HU, possibly because of what 
could be called interferences between resonant modes and 
modes in the slow manifold. 

The Coriolis force does not affect the kinetic energy 
balance nor does it modify the nonlinear part of the ex- 
act law stemming from energy conservation [T3 , [H[ when 
stated in its anisotropic version. Similarly, helicity con- 
servation is not altered by the Coriolis force, nor is (in the 
structure function formulation) the nonlinear part of the 
exact law stemming from that conservation [14| : indeed, 
helicity is conserved in the presence of rigid body rota- 
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tion and, when using structure functions as in [14j |, the 
constant rotation vector drops from the dynamical equa- 
tion for the second order moment. The fact that uniform 
rotation affects odd-order moments of the velocity but 
not even ones (at least in the linear limit of negligible 
nonlinear terms) shows that its effects are subtle. How- 
ever, it has been documented in the literature (see e.g. 
[IH) that a specific model is needed in Large Eddy Simu- 
lations (LES) of turbulent rotating flows in order to take 
into account the slowing down of energy decay because 
of waves [l6| , and the anisotropy of integral length scales 
and of dissipation. 

Furthermore, wave resonant theories are non uniform 
in scale and it is not clear whether their predictions are 
verified in high Reynolds number flows in the laboratory, 
in the environment or in direct numerical simulations 
(DNS). In the decaying case there may be differ- 
ent temporal regimes in which different mechanisms pre- 
vail ■ In a forced flow, high resolutions and long-time 
integrations are needed in order to resolve the different 
spatial regimes that may develo p (e .g. inverse and di- 
rect cascades of energy (IH [l9l. |20L l2ll |22| and direct 
cascades of helicity |23|), as well as the short-time wave 
regime versus the long-time turbulent regime. 

Modeling becomes a necessity in order to reach an un- 
derstanding of these flows at high Reynolds numbers as 
occurs in astrophysics and geophysics. Using a closure 
set of integro-differential equations for the energy spectra 
(see and references therein) or the weak turbulence 
framework developed in [J] is a powerful tool for suffi- 
ciently small Rossby number, but two difficulties have 
to be overcome. On the one hand, such techniques are 
valid in the limit Ro — > and yet the flows one tries to 
model may present inhomogeneities (in space, in time, in 
scale) that are affected differently by the rotation. On 
the other hand, the complexity of the so-called weak tur- 
bulence kinetic equations and in particular their depen- 
dence on the angle between the wavevector (in a Fourier 
decomposition) and the rotation axis (to be taken as the 
z axis in what follows) necessitates a regular discretiza- 
tion in angle as opposed to an exponential discretization 
in wavenumber, the latter working because of the self- 
similarity of the known (power law) spectral solutions to 
the equations. This angular dependency makes the clo- 
sure or weak turbulence equations difficult and costly to 
use (see however flol|). 

The anisotropy of a rotating flow comes from the non- 
linear terms (and resonant triadic interactions) [l6[, at 
least for one-time second order statistics (24|, because of 
a loss of phase information. An explicit example of an 
initial condition that gives rise to elongated vortices is 
given in (l6| ; these authors predict linear growth of the 
integral scale, as opposed to the classical Kolmogorov 
~ t 2 ' 7 law, and as observed in laboratory experiments 
. This anisotropy is also linked to an asymmetry (pre- 
dominance of cyclonic event over anti-cyclonic) for times 
of order I/O, as measured for example by the skewness of 
the vertical component of the vorticity [25|,[26|]. However, 



it has been shown by several authors that the expected 
bi-dimensionalization of the flow is only realized partially, 
and small scale eddies may not follow such a dynamics; 
in which case, one expects the small scale eddies (i.e., 
those that are to be modeled in an LES approach) to be 
somewhat isotropic. It may thus be envisageable to use, 
as a model of small scales, a methodology developed for 
isotropic flows. 

It is in this context that we extend the spectral model 
derived in [27j to the case of forced rotating flows, com- 
paring the results of the model to high resolution DNS 
[23 1 f° r forced rotating turbulence down to Ro ~ 0.03. 
The model is based on the EDQNM closure to compute 
eddy viscosity and eddy noise. It adapts dynamically to 
the inertial index of the energy spectrum, and as a result 
it is well suited to study rotating turbulence for which the 
scaling laws are not well known, and may change with the 
Rossby number, and also (at fixed Rossby number) as the 
system evolves and an inverse cascade develops. The next 
section poses the problem in terms of equations and mod- 
els and gives the numerical set-up. We then describe the 
results for the isotropic LES model, examining energetic 
balance, structures, spectra and higher-order statistics. 
Finally, the last section presents our conclusions. 

II. EQUATIONS AND SPECTRAL MODELING 

A. Primitive equations 

The dynamical equations can be written in terms of the 
Fourier coefficients of the velocity field defined as usual 
as: 

v(k, t) = JJJ°° v(x, t)e- jk x dx . (1) 

In the rotating frame, and including the centrifugal force 
in the pressure term, the equations are: 

(J^ + vk2 ) v o t {k,t)+29,P a pe l}zi u 1 {k,t) = t a (k,t)+F a {k) 

(2) 

together with the incompressibility condition k • v = 0; 
v is the kinematic viscosity, F(k) is the Fourier trans- 
form of the forcing function, P a p — 5 a f} — kakp/k 2 is the 
projection operator, f2 is the rotation rate and t(k, t) is 
a bilinear operator for the kinetic energy transfer written 
as: 

t a (k,t) = -iP a/3 (k)fc 7 ^2 W /3(P>*)«7 (<!>*) ■ ( 3 ) 

p+q=k 

Note that P a p is a projector that allows us to take the 
pressure term of the velocity equation into account via 
a Poisson formulation and ensures that the velocity re- 
mains divergence-free including in the presence of rota- 
tion. Finally note that the total energy = (v 2 /2) 
and the helicity (v ■ u) (with ui = V x v) are invariants 
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of the three-dimensional equations in the ideal case, i.e. 
in the absence of viscous dissipation [y = 0). Besides the 
Reynolds number and the Rossby number defined pre- 
viously, one can also introduce dimensionless numbers 
based on small-scales as produced by the turbulent flow; 
the simplest way to do that, traditionally, is to base such 
parameters on the vorticity through the Taylor scale A 
defined as: 



A = 2tt 



jE{k)dk 
jE{k)k 2 dk 



1/2 



(4) 



the Taylor Reynolds number is then: 

v 

One can also define a quantity called the micro-Rossby 
number pll | which is useful to determine the regime of 
the small scale turbulence and the slope of the energy 
spectrum [25|. It reads: 



ROu, = 



■in 



(5) 



where uj rms stands for the r.m.s. vorticity; note that it 
is proportional to the Rossby number evaluated at the 
Taylor scale. Finally, we also define the Ekman number: 



Ek 



Ro 
~Re 



20c 2 ,' 



(6) 



where Lq = 2-7r/fco is the scale associated with the forcing 
at kg — 2. The direct numerical simulations used in this 
paper (runs Id, lid and Hid respectively, see Table [TJ) in 
order to assess the validity of the LES are those labeled 
A3, A4 and A6 respectively in [22j (hereafter, paper I). 
For all these runs, the forcing function is a Taylor- Green 
(TG) vortex with amplitude Fq: 



F x = 

Fy = 

F z = 



F sin(fcoa;) cos(koy) cos(koz) 
—Fq cos(fco£) sin(fcoy) cos(fcoz) 
; 



(7) 



the third component of the forcing is equal to zero but 
the velocity in the z-direction grows through nonlinear 
interactions. Moreover, the forcing injects no energy in 
modes with k z = 0, and as a result any amplification 
observed in strongly rotating cases must be only due to 
a cascade process. Finally, the resulting flow has a small 
spectral anisotropy with slightly more energy in the z 
direction [22j , an effect which is the opposite of the ten- 
dency towards two-dimensionalization that develops in 
rotating turbulence. 

The numerical computations using the above forc- 
ing arc thus either Direct Numerical Simulations of the 
Navier-Stokes equations with 256 3 grid points, or Large 
Eddy Simulations on grids of 64 3 points; the axis of ro- 
tation is the z-axis, and the flow is initially at rest. Note 



that the TG flow is widely used in experimental devices 
to study turbulence and its effect on the generation of 
magnetic fields [28[ even though the TG vortex has no 
net helicity due to its symmetries; because of this latter 
property, the LES model used here will not include the 
helicity eddy viscosity derived in [27} (Paper II hereafter) . 
The turnover time at the forcing scale is then defined as 
t~nl = Lq/Uo where Uq = y/ (v 2 ) is the r.m.s. velocity 
measured in the turbulent steady state as stated previ- 
ously, at the onset of the inverse cascade at low Rossby 
number. Note that the amplitude of the forcing F in 
each simulation is increased with Q to have Uq ~ 1 in all 
the runs before the inverse cascade sets in (see [22[ for 
more details on the DNS runs). 

Finally, as the issue of the direction of the energy cas- 
cade (direct and/or inverse) is an important issue in ro- 
tating turbulence, a useful diagnostic in this context is to 
examine the behavior of the skewness (normalized third- 
order moment corresponding to energy transfer) based 
on the velocity derivative, 



S k = 



((d x v x ) 3 ) 



{{d x v x f) 



2\3/2 



(8) 



B. The isotropic EDQNM closure 

The Large Eddy Simulation model (LES) derived in 
[27j for non-rotating Navier-Stokes flows is now extended 
to the rotating case in its non-helical version (LES-P of 
Paper II). In other words, intrinsic variations of the he- 
licity spectra are not taken into account in the present 
work in the evaluation of the transport coefficients used 
in our LES model. The first step of the model is to em- 
ploy a spectral filtering of the equations; this operation 
consists in truncating all velocity components at wave 
vectors k such that |k| = k > k c , where k c is a so-called 
cutoff wave number. Since the scales associated with k c 
are presumably much larger than the actual dissipative 
small scales in a high Reynolds number flow, one needs 
to model the transfer between the large (resolved) scales 
and the small (subgrid unresolved) scales of the flow. In 
order to approximate these transfer terms, the behav- 
ior of the energy spectrum after the cutoff wave number 
has to be estimated. We therefore define an intermediate 
range, lying between k c and 3fc c , where the energy spec- 
trum is assumed to present a power-law behavior possibly 
followed by an exponential decrease: 



E v {k,t) 



-5%k 



(9) 



The coefficients a^, and Eq are computed at each 
time step, through a mean square fit of the resolved en- 
ergy spectrum. In a second step, one can write the fol- 
lowing model equations (omitting forcing): 



[d t + (u(k\k c ,t) + u) fc 2 ]« Q (k,i) 
= t<(k I t)-2nP a/J e / j^u y (k > t) 



(10) 



4 



where the < symbol indicates that the nonlinear trans- 
fer terms are integrated over a truncated domain defined 
such that p + q = k with |p| = p, |q| = q < k c . The eddy 
viscosity v (k\k c , t) is expressed as 



u(k\k c ,t)- 



A> 



SE 2 (k,p,q,t) 
• 2k 2 E v (k 1 t) 



dpdq 



The function Se 2 (k,p, q, t) corresponds to the so-called 
absorption term (linear in the energy spectrum E v (k, t)) 
in the EDQNM nonlinear transfer, lending itself in par- 
ticular to an expression for the turbulent eddy viscosity, 
as is well known; A > is the integration domain over (k, 
p, q) triangles, such that p and/or q are larger than k Cl 
and both p and q are smaller than 3fc c . 

Finally, to take into account the effect of the emis- 
sion (eddy-noise) term involved in the EDQNM nonlin- 
ear transfer (Se 1 (k,p, q, t)), we use a reconstruction field 
procedure which enables us to partly rebuild the phase 
relationships between the three spectral components of 
the velocity field, as explained in detail in Paper II [27| . 
The functions S , B 1 (fc,p,(7,t) and SE 2 (k,p,q,t) appearing 
in the transport coefficients used in the LES are writ- 
ten for completion in the Appendix. Note that, although 
isotropic, the subgrid model we use in this paper has an 
important feature: it adjusts dynamically to the energy 
spectrum instead of assuming a given spectral law, usu- 
ally the classical Kolmogorov law, E{k) ~ fc~ 5 / 3 . This 
allows for exploration of flows for which a theory to pre- 
dict spectral indices is not available. Also note that the 
reconstruction procedure differs as well from traditional 
implementations insofar as it tries to keep the phase in- 
formation of the small-scales. 



III. ROTATION AND ISOTROPY 

One of the effects of rotation on a flow is to induce 
anisotropy, as in the formation of large-scale columnar 
vortices. In that light, we explore in this section the 
anisotropic properties of a DNS at low Rossby number to 
see whether or not it is relevant to use a model based on 
isotropic assumptions to simulate a flow subjected to ro- 
tation. The LES model we propose to use approximates, 
as is customary, the transfer from the large to the small 
scales, but most of the modeled interactions are between 
small scales because of the value of k c (chosen to be in 
all cases larger than the energy injection wavenumber), 
and because most of the modes in a turbulent flow are 
in the small scales (recall that the number of modes in a 
given isotropic shell fc, varies as kf). 

We therefore investigate now the properties of the 
small scales of flows forced with the Taylor- Green vor- 
tex (see Eq. [7|) and subjected to rotation, with fco = 2 
and at a Rossby number Ro = 0.03; we perform a DNS 
on a grid of 256 3 points and with the flow being initially 
at rest. To measure anisotropy, we introduce two differ- 
ent quantities, a spatial one and a spectral one, denoted 
respectively I D (for dimensional) and I c (for Craya [2!| ; 



see also [301). Another measure of anisotropy linked to 
the so-called polarization anisotropy, following [3lJ (see 
also [IH), is discussed later in Section ITVB I The spatial 
coefficient I D evaluates the averaged ratio between the 
intensity of the velocity in the perpendicular direction 
V±(x,t) and in the parallel direction Vj|(x, i), with J_, || 
referring to the z-axis of rotation. The velocity field can 
be expressed as a function of these two components as 
v(x, t) = Vj|(x, f)e|| + Vj_(x, t), where e|| is the unit vec- 
tor associated to the axis of rotation and Vj_(x, t) is the 
velocity field projected on the plane perpendicular to ey. 
The spatial anisotropy coefficient therefore reads: 



tD 



V||(x,i) 



(11) 



The spectral coefficient I c is computed as in (34|: for 
each wavevector k, an orthonormal reference frame is 
defined as (k/|k|, ei(k)/|ei(k)|, e 2 (k)/|e 2 (k)|), with 
e^k) = k x z and e 2 (k) = k x ei(k), where z is the 
vertical unit wavevector. In that frame, since the in- 
compressibility condition yields k- v(k) = 0, v(k) is only 
determined by its two components Vi(k) and v 2 (k). This 
second anisotropy coefficient is then defined as 



y/(\vi(W)/(W2(V\ 2 ) 



(12) 



Both I D and I c are such that they have unit values for 
fully isotropic flows. 

In Fig. [T]we plot the temporal evolution of the total en- 
ergy, the time being expressed in units of the eddy turn- 
over time. Note the long interval before turbulence fully 
develops, as rotation is strong and the run was started 
from a fluid at rest. Indeed, before the energy starts to 
grow at t w 90, one can observe a long transient during 
which the energy displays damped oscillations in time 
(see Paper I). This transient is linked to the effect of 
rotation and its duration increases linearly with f2, i.e. 
as the inverse of the Rossby number. During this first 
stage, the energy dissipation rate is small and the energy 
spectrum is very steep. Later, at f w 90, the enstrophy 
starts to grow and the energy dissipation rate increases. 
The energy also grows and an inverse cascade of energy 
develops. Turbulence sets in and the small-scale energy 
spectrum develops an inertial range with scaling close to 
~ fcj 2 (see Paper I for more details). 

In order to quantify the importance of anisotropy at 
what would be the sub-grid scales in a LES of rotating 
flows, we start by noting that the velocity (in particu- 
lar when an inverse cascade of energy develops at small 
enough Rossby number) is dominated by the large scales 
whereas the modeling will occur in the small scales of 
the flow. In this context, we introduce a band-pass fil- 
ter of the DNS data in order to concentrate the analysis 
on small-scale properties of the flow. The filtered field is 
given in Fourier space by all the velocity components at 
a wavector k such that 32 < |k| < 64; note that, for this 
DNS using a classical 2/3 dealiasing rule, the maximum 
wavenumber is k max — 85. As a result, the band-pass 
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FIG. 1: Temporal evolution of the energy for the flow in a 
DNS with Ro = 0.03; it displays two temporal phases, one 
dominated by wave interactions and one, for t larger than 
100, corresponding to an inverse cascade of energy. 



filter can be interpreted as preserving the small scales of 
the direct cascade inertial range. 

Figure [3] represents the time history of the I c and I D 
anisotropy coefficients for the complete DNS (dash line) 
and for the band-pass filtered velocity fields of the flow 
at Ro = 0.03 (ovals). The Craya spectral coefficient I c 
of the complete DNS field remains close to unity during 
the whole simulation, indicating that globally the flow 
is close to an isotropic state. However, the directional 
coefficient I D exhibits three different regimes in the full 
DNS: a first phase between t = and t ~ 40 during 
which it oscillates, with an amplitude that decreases with 
time, and a second phase, between t ~ 40 and t ~ 90, 
with this coefficient remaining constant at a value close to 
unity, meaning that no direction is privileged in the flow. 
Finally, in a third and last phase, which begins when the 
turbulence starts to develop, I D strongly increases with 
time. This behavior is the signature of the generation 
of intense columnar structures within the flow, within 
which the perpendicular component of the velocity field 
dominates the parallel one. 

The behavior of these coefficients is completely differ- 
ent for the filtered, small-scale, field. Indeed, the small 
scales are strongly anisotropic during the transient pe- 
riod before the turbulence develops, with a maximum 
value for I c of the order of 3 (and 5 for I D ). In this 
phase, the directional anisotropy coefficient clearly shows 
that the perpendicular component of the velocity domi- 
nates the parallel one, and therefore that the small scales 
are mostly bi-dimensional. At t ~ 80, both coefficients 
drop rather abruptly to a value of order unity, indicat- 
ing that when the turbulence develops the small scales 
become isotropic corresponding to a standard cascade of 
energy to small scales (note that the scales for which 
the anisotropic and inverse accumulation of energy takes 
place are eliminated by the band-pass filter). 

With this study of the small-scale behavior of a flow 
subjected to moderate rotation, we see that an isotropic 
LES model cannot be used to treat every phase of the 
flow. Indeed, in the early transient phase, a model based 
on isotropic assumptions will not be able to approxi- 
mate properly the transfer between the subgrid scales 
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FIG. 2: Temporal evolution for a flow with Ro = 0.03 of 
the Craya anisotropy coefficient I c (top) and the directional 
anisotropy coefficient I D (bottom) (see eqs. 111! I12[) . for 
the full DNS velocity field (dash line) and the filtered DNS 
field (ovals) defined as a band-pass filter for wavevectors in 
\K\ 6 [32,64]. Note the sharp transition towards isotropy of 
the small scales for t ~ 100, as both the direct and inverse 
turbulence cascades finally develop. 



and the resolved scales. We therefore decide to only 
use our model to study the turbulent regime of rotating 
flows, after t « 100 in the case of Figs. [TJ [21 Moreover, 
this is consistent with the fact that a LES is designed 
to study turbulent flows, and cannot handle transitional 
(laminar, wave-dominated) flows. In the case of rotating 
flows starting from a fluid at rest, turbulence only devel- 
ops after a transient time that depends linearly with the 
magnitude of the rotation. Note that in many studies, 
simulations of rotating flows are started from a previous 
turbulent steady state, and in that case our LES should 
have no problem to adapt as the spectral index changes 
with the evolution of the system. 

Note that both coefficients I c and I D are relevant 
quantities in the context of this EDQNM-based LES: 
the behavior of I c justifies the assumption of "spectral 
isotropy" (i.e., dealing with k instead of (ku, k±) at small 
scales); on the other hand, the behavior of I D justifies 
the isotropic reconstruction done with the eddy-noise, 
because I D « 1 is a measure of variance isotropy. 



IV. NUMERICAL TESTS OF THE LES 

We now test our LES model against direct numerical 
simulations with different Rossby numbers. As stated 
before, the forcing used is the Taylor-Green vortex (see 
Eq. [7]) at fco = 2. For each simulation, we follow the nu- 
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merical procedure described in Paper I; namely, we vary 
the rotation rate leading to three different Rossby num- 
bers: 0.03, 0.17, and 0.35. The simulation parameters 
are summarized in Table HI The flow evolves in a peri- 
odic box, with 256 3 grid points for the DNS and 64 3 grid 
points for the LES. The "reduced-DNS" results, in the 
table and figures, refer to the filtered DNS data on a grid 
of 64 3 grid-points, corresponding to the limited informa- 
tion contained in the LES grid. Since we are interested in 
studying only the modeling of the turbulent regime, we 
start the LES simulations from the reduced-DNS data at 
a time after the end of the transient phase. However, 
if the LES is started from a fluid at rest (i.e., started 
like the DNS at t = 0), no significant differences are ob- 
served with the procedure of starting the LES at the end 
of the transient phase, except that the transient regime 
in the flow with Ro = 0.03 is shorter. This accelerated 
evolution of the LES at low Rossby number during the 
transient when compared to the DNS can be easily ex- 
plained considering the inclusion of transport coefficients 
in the LES which assumes that a turbulent flow is already 
present. 



A. Global behavior of the flow 

The first test of the model is to examine the tempo- 
ral evolution of the flow. This is displayed in Fig. [3] 
for the three Rossby numbers analyzed. The overall be- 
haviors of the DNS and of the LES are similar in ampli- 
tude and in time scales. At intermediate Rossby num- 
bers {Ro = 0.17), the precise evolution of the DNS is 
not followed although the energy obtained with the LES 
remains close to the DNS one. For the simulation at 
Ro = 0.03, an inverse cascade develops after t ~ f20 
leading to a strong increase of the total energy. Although 
the LES model does not take wave interactions explicitly 
into account, it allows to reproduce this transfer of en- 
ergy from the small scales to the very large ones with 
good accuracy; indeed, a scaling argument shows that in 
the small scales, the eddy turn-over time is shorter than 
the time associated with waves and nonlinearities prevail. 
The LES is taking into account the interactions with the 
waves in an implicit way by changing the EDQNM time 
scale dynamically with the slope of the energy spectrum 
at large scales; this could be interpreted as "reversed" 
Kraichnan-likc phenomenology. Note that the run at in- 
termediate Rossby number has higher values of the en- 
ergy because the forcing amplitude is larger than for the 
run at Ro = 0.35. 

When looking at the time-averaged isotropic energy 
spectra (see Fig. HJ for the two flows at the largest Rossby 
numbers, one can see that a good agreement is obtained. 
This figure also allows us to better understand the differ- 
ence in the temporal evolution of the energy computed 
from the DNS and the LES data at Ro = 0.17 (see Fig. 
|3]). Indeed, although the model gives a good estimation 
of the DNS spectra at small scales, at very large scale 
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FIG. 3: Temporal evolution of the total energy for runs Id 
(DNS: 256 3 ) and IL (LES: 64 3 ) at Ro = 0.03 (top), runs lid 
(DNS: 256 3 ) and IIL (LES: 64 3 ) at Ro = 0.17 (middle) , and 
runs Hid (DNS: 256 3 ) and IIIL (LES: 64 3 ) at Ro = 0.35 
(bottom). DNS runs with dash line, LES runs with solid line. 
Note the change of values on both axes for the low Rossby 
runs (top) because of the delay in the development of the 
turbulent phase, when the LES is started, and the ensuing 
accumulation of energy due to the inverse cascade now taking 
place at that low Rossby number. 



(and particularly at k = 2) non-negligible differences ap- 
pear with the DNS, differences to which the total energy 
is sensitive. Note that a smaller difference between LES 
and DNS spectra can be observed at k = 2 for the run 
at the higher Rossby number, Ro = 0.35. Otherwise, 
the spectrum is well approximated by the LES at all the 
other scales. 

Similarly, when decomposing the energy spectra into 
their perpendicular and parallel components, a good 
agreement is reached at large scales, except again at 
k = 2 for the perpendicular spectrum of the flow at 
Ro = 0.17 (see Fig. [5]). On the contrary, at small scales, 
the model seems to underestimate the spectra obtained 
by the DNS. This behavior is in fact due to the differ- 
ence in resolution between the DNS and the LES: as k± 
and ku increase, the difference between the amount of 
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TABLE I: Parameters of the simulations: Resolution N 3 , Rossby number Ro based on the forcing scale Lo, Taylor microscale A 

and integral scale L, r.m.s.velocity Uo = (v 2 ) 1 '' 2 , integral Reynolds number Re = UqL/v and eddy turnover time tnl = Lo/Uo; 
t m is the final time of the computation. Note that the r label in the nomenclature of the runs stands for reduced data obtained 
by filtering in spectral space to 64 3 points the original 256 3 DNS data A, L, Re and tnl are evaluated at the final time of the 
simulation for runs I which undergoes an inverse cascade, whereas they are averaged during the stationary phase of simulations 
II and III which are at higher Rossby numbers and do not undergo any significant inverse energy transfer. 
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FIG. 4: Time averaged energy spectra for runs lid (DNS: 
256 3 ) and IIL (LES: 64 3 ) at Ro = 0.17 (top), and runs Hid 
(DNS: 256 3 ) and IIIL (LES: 64 3 ) at Ro = 0.35 (bottom). 
Time averaging is performed from t = 25 to t = 45 for both 
simulations. Note the good agreement except possibly near 
k = 2 corresponding to the forcing scale, indicative of a lack 
of adjustment of the LES at that scale, in particular for the 
perpendicular spectra, see Fig. [5] 



modes taken into account in the evolution of these spec- 
tra for the DNS and for the LES increases as well. Note 
that the fc|| shells have the same number of modes in- 
dependently of the value of kii (they are planes), while 
the number of modes in the k± shells grow as k± (they 
are cylinders), and this number grows as k ' 1 in dimen- 
sion D for isotropic (spherical) shells. We have checked 
that, when making the comparison between the LES and 
the reduced DNS for instantaneous spectra, the discrep- 





DNS 

+ LES 




FIG. 5: Time averaged parallel (left) and perpendicular 
(right) energy spectra, for runs lid (DNS: 256 3 ) and IIL 
(LES: 64 3 ) at Ro = 0.17 (top), and runs Hid (DNS: 256 3 ) 
and IIIL (LES: 64 3 ) at Ro = 0.35 (bottom). 



ancy observed at high wavenumber disappears. 

As mentioned earlier, the micro-Rossby number mea- 
sures how strong the imposed rotation is in the flow at 
the Taylor microscale, when compared to the r.m.s. vor- 
ticity developed by the turbulence. Its time evolution is 
shown in Fig. [5] for all runs. Because the micro-Rossby 
number emphasizes small scales that are not all present in 
an LES, i?o w is also computed in the reduced-DNS. We 
observe a good agreement between the truncated DNS 
and the LES, although the model slightly underestimates 
Ro u for the two simulations at larger Rossby number. 
This behavior can be explained by enstrophy production 
in the LES, and the backscattering of energy from sub- 
grid scales to resolved scale associated to the eddy noise, 
which is perhaps not strong enough. 

In Table [II] we give the values of the characteristic 
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FIG. 6: Temporal evolution of the micro Rossby number 
Rou: in flows with Ro = 0.03 (top), Ro = 0.17 (middle) and 
Ro = 0.35 (bottom), for DNS on 256 3 grid points (dash line), 
filtered data of the DNS to a 64 3 grid (triangles) and LES 
(solid line) on a 64 3 grid. Note again the different scale on 
the axes, in particular for the lowest Rossby number in which 
case Rouj approaches unity as the inverse cascade develops 
and energy and turbulence intensity grow. 



parallel and perpendicular integral length scales (respec- 



tively L» and L±) denned as: 



in = 



Lj 



E(k {l )k^ 1 dk ll 



It 



ax E(k\\)dk\\ 
E(k 1 _)kJ 1 dk± 



E(k ± )dk ± 



(13) 
(14) 



and computed at the final simulation time of each flow 
(note that the fc|| = mode is not included in the defi- 
nition) . Even if the values obtained by the LES data do 
not exactly correspond to the DNS values, they remain 
close; their difference can be explained by the same ar- 
gument evoked before on the slight discrepancy between 
LES and DNS parallel and perpendicular energy spec- 
tra. Note that the perpendicular length scale is signifi- 
cantly larger for the lowest Rossby number, but the par- 



TABLE II: Characteristic integral length scales L± and Ly 
measured at different times t m for the three different Rossby 
numbers studied in this paper. Note that, at the lowest 
Rossby number (Ro = 0.03, runs I), the perpendicular inte- 
gral length scale is significantly larger than for more moderate 
rotation, because of the inverse cascade. 
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allel length scales are comparable in all three runs. This 
is linked to the fact that the inverse cascade of energy 
which takes place at low Rossby number is dominated by 
quasi-two-dimensional modes; the parallel spectrum does 
not undergo an inverse cascade, although energy docs pile 
up at fcii = mainly through resonant coupling of waves. 



B. Measures of anisotropy 

Rotating flows are known to develop anisotropics and 
we now turn our attention to this point. In order to esti- 
mate the anisotropy of the different flows, we use the co- 
efficients I D and I c defined earlier in Eqs. (fTTj) and (fl"2j) . 
They are shown as a function of time in Fig. [7] for the 
DNS (dash line), the reduced DNS data truncated to the 
LES resolution (triangles), and the LES (solid line) with 
Ro = 0.03. A very good match can be observed between 
the Craya coefficient I c computed from the rcduccd-DNS 
data and the one computed with the data from the LES 
model, whereas the coefficient computed with the full 
DNS data evolves on a lower level than the two other 
ones. This is due to the fact that the small scales of 
the field (i.e. scales with k > k c ) are taken into account 
in the spatial averaging process we perform to calculate 
this coefficient. We saw in Section Mil that these small 
scales arc more isotropic with a corresponding coefficient 
I c near unity, so when they are taken into account in 
the computation of the Craya coefficient they lower its 
value. The small scales in the DNS are more isotropic, 
and as a result, the LES flow, which preserves a smaller 
amount of these scales, is globally more anisotropic and 
has a larger value of this coefficient. 

As already observed in Fig. [51 the directional coef- 
ficient I D is strongly dominated by the large scales of 
the field, such as columnar structures appearing in the 
flow as a result of the inverse cascade process. Therefore, 
when we compare the time history of this coefficient for 
the DNS and the reduced-DNS, no noticeable difference 
appears. Once again our LES model predicts very well 
the evolution of this coefficient, even though the perpen- 
dicular component of the velocity clearly dominates over 
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FIG. 7: Evolution of the spectral and directional isotropy 
coefficients I c and I D (see text) for runs Ir (reduced DNS: 
64 3 ) and IL (LES: 64 3 ) at low Rossby number (Ro = 0.03). 
Information about the lower rotation cases is given in Table 
IIHI Isotropy obtains when both coefficients are close to unity, 
and we note that the directional coefficient, related to real- 
space structures, indicates a stronger departure from isotropy 
than when measuring anisotropy in Fourier space as Ic does 
(see Eqs. 1 1 1 1 I12|l , once the inverse cascade sets up and strong 
columnar vortices develop. Larger I D also obtain for the runs 
at lower Rossby numbers (see Table tTTT|) . 



TABLE III: Craya and directional isotropic coefficients I 
and I D for the simulations at Ro = 0.17 and Ro = 0.35. 
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the parallel one. We also note that the model allows 
for a good estimation of both these coefficients for the 
simulations at larger Rossby numbers (Ro = 0.17 and 
Ro = 0.35), as shown in Table UTTl 

In our investigation of anisotropy of rotating flows, we 
finally study the behavior of the bij anisotropy tensor de- 
fined below (see e.g. [HI for reference); it is linked to the 
so-called "polarization" anisotro py i ntroduced in [3l| and 
as also discussed in [32[ (see also j33|). This tensor, which 



FIG. 8: Temporal evolution of the b zz component of the polar- 
ization anisotropy tensor (see Eq. I15p for runs with Ro = 0.03 
(top) and Ro — 0.35 (bottom), in the former case at late times 
once the pile-up of energy at large scales has begun. The DNS 
are plotted with a dash line, the reduced DNS, truncated to 
the resolution of the LES, are given with triangles and the 
solid line corresponds to the LES. Note again the tendency to 
a two-dimensional state at late times, with b zz — > —1/3, for 
the low Rossby number runs. 



is defined as: 



b - ^ 
R>u 



l 3 

3 ' 



(15) 



is based on the Reynolds stress tensor Ri 



) u i( x )), 



with summation upon the subscript /. In Fig. [5] we rep- 
resent the temporal evolution of the b zz component of 
the anisotropy tensor for runs I and III, at respectively 
Ro = 0.03 and Ro = 0.35. We first notice that the LES 
model predicts well the evolution of this coefficient for 
both simulations. Secondly, the development of a pre- 
ferred direction in the flow at Ro = 0.03 (already ob- 
served in Fig. [7] through the increase of the directional 
coefficient I D in the inverse cascade), is also visible in this 
figure. Indeed, b zz tends to —1/3 as time increases, since 
u z (x) becomes negligible when compared to the horizon- 
tal components v x (x) and v y (x) . 



C. Statistical analysis 

In this section, we investigate the statistics of the sim- 
ulations at Ro = 0.03. Instantaneous probability density 
functions (or PDFs) of the longitudinal and transverse 
derivative of the x-component of the velocity field are 
computed and plotted in Fig. [!|]at time t — 132, in the 
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inverse cascade. The PDFs computed on the full DNS 
data, the reduced-DNS, and the LES, agree well for the 
case of the longitudinal derivative. In the case of the 
transverse derivative, the DNS data differ from both the 
LES and the reduced-DNS data, the latter two display- 
ing wider wings and being almost superimposed. It is 
well known that the small scales of a flow may have a 
strong influence on the distribution of velocity deriva- 
tives with strong velocity gradients appearing at small 
scale, and that transverse derivatives show stronger tails 
in the pdfs (and therefore enhanced intermittency) than 
longitudinal derivatives. It is not clear whether this is 
the effect of more sensitivity to the intermittency in the 
transverse increments or whether it is the effect of small- 
scale anisotropy, but since the differences are stronger for 
the velocity derivatives taken in the direction of rotation, 
it may be attributed to anisotropics. 




locity derivative, i.e. its normalized third order moment. 
The skewness, which measures the departure from Gaus- 
sian statistics, is usually negative for the longitudinal 
derivatives of a turbulent flow and oscillates around zero 
for the lateral ones. In Fig. [10] we show the time his- 
tory of iSfc. As for the energy, the LES model gives a 
correct prediction of the skewness for 86 < t < 145, al- 
though around t ~ 140, some discrepancy can be found 
that could be associated with the development of struc- 
tures. Note that this difference can be also associated 
with a slight discrepancy in total energy at around the 
same time (see Fig. [3]). 
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FIG. 10: Temporal evolution of the skewness for the longi- 
tudinal velocity derivative dv x /dx (see Eq. [8}, for runs Ir 
(reduced DNS: 64 3 ) and IL (LES: 64 3 ) at Ro = 0.03. 
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FIG. 9: Probability density function of the velocity deriva- 
tives dvx/dx (top), dv x /dy (middle), and dv x /dz (bottom), 
for runs Id (DNS: 256 3 ), Ir (reduced DNS: 64 3 ), and IL (LES: 
64 3 ) at Ro = 0.03 and t = 132. As usual, dash line is for the 
full DNS flow, triangles for the reduced (truncated) DNS and 
solid line for the LES. 

In order to quantify the distributions of velocity fluc- 
tuations and the differences between DNS and LES data, 
we now compute the skewness Sk of the longitudinal ve- 



D. Visualization in physical space 

We finally present a visualization in physical space of 
the velocity intensity at t = 132 for the flow at Ro = 0.03. 
At this time of the simulation the inverse cascade already 
took place and most of the flow energy was transferred 
to the fcii = plane. We noted earlier that the TG 
flow injects no energy in the k± = 1 shell nor in the 
fen = shell. So all energy we see at large scale is the 
result of inverse cascade (in the former case) and of two- 
dimensionalization (in the latter case). The evidence for 
the inverse cascade in this paper is given by the time evo- 
lution of the energy in Figs. 1 and 3 (see also Paper I, 
where fluxes arc studied in detail). The accumulation of 
energy in this plane leads to the formation of columns as 
can be observed in Fig. [TT] Although the structures are 
quasi-bidimensional, the isotropic LES model allows to 
reproduce them quite correctly. The spatial position dif- 
fers slightly from the structure obtained by the DNS, but 
its size and intensity are well approximated. When exam- 
ining the temporal evolution of the maximum of velocity 
(not shown), a good agreement occurs at all times. Note 
that this is a forced run visualized after « 130 turnover 
times; as a result of the intrinsic sensitivity of turbulent 
flows due to their inherent unpredictability after a Lya- 
pounov time of the order of a few turn-over times, the 
spatial position of the structures is not expected to be 
reproduced exactly by the LES. 
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FIG. 11: Volume rendering of the velocity modulus for flows 
at Ro = 0.03 - runs Ir (Reduced DNS: 64 3 , top), and IL 
(LES: 64 3 bottom) - at time t — 132 for both simulations. 
The flow is dominated by a large eddy, but smaller vertical 
structures are visible as well. 



V. CONCLUSION 

We present in this paper a Large Eddy Simulation 
model for high Reynolds number rotating flows using 
a previously derived sub-grid model [27| (see also (3^ |) 
based on the isotropic EDQNM two-point closure with 
eddy viscosity and eddy noise. We show that, down 
to Rossby numbers of 0.03, the small-scales are suffi- 
ciently isotropic for the model to perform reasonably well. 
There are numerous laboratory experiments with which 
the comparison presented here could be extended, follow- 
ing the work in [2l[ using several EDQNM-based closures 
(see also [37j for the stratified anisotropic case) . The ad- 
vantage of using two-point closures such as EDQNM as 
a model for turbulent flows in the presence of rotation is 
that it allows for building scaling laws at a relatively low 



computational cost and with the possibility of doing ana- 
lytical estimations of nonlinear transfer (sec for example 
[2l|). The model presented in this paper is much simpler 
since it is built on the isotropic three-dimensional ver- 
sion of the EDQNM; it is thus more limited in its scope 
insofar as it may not be able to explore very low Rossby 
numbers. On the other hand, following the standard LES 
methodology with spatially resolved large scales, and tur- 
bulent coefficients to model the sub-grid fluctuations, it 
allows to access more detailed features of the flows such 
as high-order statistics (e.g., PDFs) as well as spatial 
structures. 

Also, the LES used in this paper adapts dynamically 
depending on the spectral index of the energy at super- 
filter (resolved) scales, and the value of the turbulent 
transport coefficients vary as a result. This is important 
in the context of rotating turbulent flows because the 
power law followed by the energy spectrum in this case 
is not quite ascertained yet and does vary with time. 
Phcnomenological and theoretical predictions of this in- 
dex, as well as several recent experiments, were reviewed 
in [25(, with experimental and numerical evidence not 
quite able yet to sort out the different models or to fully 
describe the parameter space (e.g., as function of the ro- 
tation rate f2, the Reynolds number, etc.). An adequate 
LES model that can adjust to the resolved energy spec- 
trum can help in this matter but more development and 
tests are needed. A reminiscent situation is found in 
magnetohydrodynamics (MHD) when coupling the fluid 
to a magnetic field in the non-relativistic limit; the to- 
tal energy spectrum obtained analytically from the weak 
turbulence limit Hi! has been observed in the magne- 
tosphere of Jupiter [40] and in DNS [!l| , but the strong 
turbulence spectrum (or spectra in case there are differ- 
ent regions in parameter space) is a matter of debate. 

Only one specific (non-helical) forcing was explored in 
the DNS-LES comparisons studied in this paper. Fur- 
ther tests are required, considering other (non-helical) 
forcing functions, as well as forcing functions that intro- 
duce both energy and helicity in the flow. In this latter 
case, the implementation of the LES as described here 
may prove insufficient and one should also consider tak- 
ing into account the spectral properties and turbulent 
transport coefficients that include the effect of helicity, 
as done in the non-rotating case in Paper II. Such an im- 
plementation can also be of interest for non- helical flows, 
because even though helicity is not a positive definite 
quantity, local helical fluctuations develop rapidly in a 
flow through alignment of vorticity and pressure gradi- 
ents [42| ■ The properties of the model in the helical case 
in the presence of rotation will be dealt with in a forth- 
coming paper. The freely-decaying case (see [HI HH [26| 
for a global perspective) needs to be examined as well 
and is left for future work. 

As a final remark, we want to stress the importance 
of developing adequate modeling of rotating (and strat- 
ified) flows, as encountered for example in the Earth 
atmosphere. It was shown recently 43j that the max- 
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imum intensity of a hurricane depends crucially on the 
(assumed) horizontal mixing length; this implies that an 
adequate treatment of the turbulence is essential in pre- 
dicting various properties of hurricanes such as its inten- 
sity or landfall localization. A run with resolution down 
to 62 meters shows strong local winds that were unre- 
solved in previous studies 44|. If the work presented 



here (as well as most of its predecessors) is far from real- 
ity for hurricane dynamical modeling (because of its lack 
of proper boundary conditions, of stratification, of mois- 
ture, ...), it nevertheless represents a first step towards 
the goal of a better understanding of geophysical flows, 
the issue here being that sufficiently high Reynolds num- 
ber, i.e. sufficient multi-scale interactions and two-way 
coupling between the small scales and the large scales 
in turbulent fluids supporting inertial (and/or gravity) 
waves, is a desired ingredient for testing LES approaches 
to geophysical turbulence. 



T(M) 



A 



(t)S E (k,p,q,t)dpdq 



(A2) 



Here A is the integration domain with p and q such that 
(k, p, q) form a triangle, and k (t) is the relaxation time 
of the triple velocity correlations. As usual [45|, kpq (t) 
is defined as : 



1 



(A3) 



where fik expresses the rate at which the triple correla- 
tions evolve, i.e. under viscous dissipation and nonlinear 
shear. It can be written as: 



fX k = vk 2 + X K ( / q 2 E(q, t)dq} 



1/2 



(A4) 
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APPENDIX A: CLOSURE EXPRESSIONS OF 
TRANSFER TERMS 



Note that Xk is the only free parameter of the problem, 
taken equal to 0.36 to recover the Kolmogorov constant 
Ck = 1-4 for a k~ 5 / 3 classical energy spectrum. The ex- 
pressions of S E (k,p,q,t) can be further explicited (with 
the time dependency of energy spectra omitted here) as: 



S E (k,p,q,t) = -b[k 2 E{q)E(p)-p 2 E{q)E{k)] 

pq 



S El {k,p,q,t) + S E2 (k,p,q,t) 



For completeness, we recall here the expression of the 
EDQNM closure equation for the kinetic energy spec- 
trum E(k, t) without helicity (note that the Coriolis term 
vanishes in the energy equation). 



{d t + 2uk 2 )E{k,t) =T(k,t) 



(Al) 



where the nonlinear transfer terms T(k,t) is expressed 
as: 



Here, S El (k,p, <?,£), and S E2 {k,p,q,t), are respectively 
used to denote the two terms of the extensive expression 
of S E (k,p,q,t). The geometric coefficient b(k,p,q) (in 
short, b in the previous expression) is defined as: 



b=?-(xy + z 3 ) , 



(A5) 



where here, x, y, z are the cosines of the inner angles 
opposite to k, p, q. 
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